#Read in data
hcpd_data <- readr::read_csv("hcpd_n762_myelin_cog_data.csv")
#Column names of myelin measurements from 7 networks
seven_nets <- c('mean_wholebrain_T1wT2w', paste0('Yeo', 1:7, '_myelin')) #the first element is whole-brain

model_list <- lapply(1:length(seven_nets), function(i){
  right_side <- '~ s(Age) + Sex + Scanner + numNavs_sum'
  left_side <- seven_nets[[i]]
  
  mf_form <- as.formula(paste0(left_side, gsub('s\\((.*)\\)', '\\1', right_side)))
  form <- as.formula(paste0(left_side, right_side))
  
  d <- model.frame(formula = mf_form, data = hcpd_data)
  
  model_name <- ifelse(i == 1, 'fit_brm', paste0('fit_brm_', left_side))
  if(!file.exists(paste0(model_name, '.rds'))){
    fit_brm <- NULL
  } else {
    fit_brm <- readRDS(paste0(model_name, '.rds'))
  }
  return(fit_brm)
})
posterior_samples_list <- lapply(model_list, get_smooth_posterior, res, eps)
results_summaries_list <- lapply(posterior_samples_list, get_smooth_summaries)

Plots for each network

for (i in 1:length(seven_nets)){
  include_points <- TRUE #overlay points on smooth plots
  
  cat(paste0('\n\n## ', seven_nets[[i]], '{.tabset}\n\n'))  
  
  cat(paste0('\n\n### Derivative\n\n'))
  cat('Age where the median posterior derivative is at its max is marked.\n\n')
  
  aderivative <- results_summaries_list[[i]]$derivative
  anageatmax <- results_summaries_list[[i]]$max_deriv_age_posterior
  asmooth <- results_summaries_list[[i]]$smooth
  age_at_max_med_deriv <- aderivative[deriv1 == max(deriv1), Age]
  
  #In order to put the smooth on the same scale as the data, we need to add back
  #in the mean of the response variable that brms uses to compute the
  #conditional effect of the smooth. The easiest way to do this is to use the
  #output of the conditional effects output. We don't really care about the effect across
  #age so we can get the output at an arbitrary age:
  
  ce <- conditional_effects(model_list[[i]], effects = 'Age', 
                            int_conditions = list('Age' = c(8)))
  #get the model frame too, so we can plot the points to confirm this works.
  mf <- brms:::model.frame.brmsfit(model_list[[i]])
  
  #We want to reference the mean of the response variable, which I've set up to
  #be in the list `seven_nets`. I'm naming the variable `response_mf_mean` to
  #indicate that this is the mean of the response as determined by brms from the
  #model frame (`mf`) that it uses to estimate the model.
  response_mf_mean <- ce$Age[[seven_nets[[i]]]]
  
  #We now need to add the response variable mean to the smooth estimates, but we
  #only need to do that for the plot of the smooth---the derivative is already
  #in the correct units (change in mylenation over age) and does not need to be
  #adjusted. See below....
  print(ggplot(aderivative, aes(x = Age, y = deriv1)) + 
    geom_ribbon(aes(ymin = lower, ymax = upper), alpha = 1, fill = apal[2]) +  
    geom_segment(x = age_at_max_med_deriv, 
                 xend = age_at_max_med_deriv, 
                 y = aderivative[Age == age_at_max_med_deriv, lower],
                 yend = aderivative[Age == age_at_max_med_deriv, upper], color = apal[[1]]) +
    geom_line() + 
    geom_line(aes(color = compared_to_max, group = 1), size = 2, alpha = .8) + 
    geom_point(x = age_at_max_med_deriv, y = aderivative[Age == age_at_max_med_deriv, deriv1], color = apal[[1]], size = 2) +
    scale_color_manual(breaks = c('less_steep', 'as_steep'), values = apal[c(4,5)], 
                       labels = c('Less steep', 'As steep'),
                       name = 'Region, as compared \nto maximum derivative, \nis...') +
    jftheme)

    cat(paste0('\n\n### Age at maximum derivative\n\n'))
    
    bayesplot::color_scheme_set('red')
    print(bayesplot::mcmc_areas(data.frame(max_age = anageatmax$max_age, Chain = rep(1:4, each = 2500)),
                          prob = 0.95,
                          prob_outer = 0.99) + jftheme)
    cat(paste0('\n\n### Difference from derivative at steepest age\n\n'))
    
    print(ggplot(aderivative, aes(x = Age, y = diff_from_max)) + 
      geom_ribbon(aes(ymin = diff_lower, ymax = diff_upper), alpha = 1, fill = apal[2]) +  
      geom_line(color = apal[[1]]) + 
      geom_line(aes(color = compared_to_max, group = 1), size = 2, alpha = .8) + 
      scale_color_manual(breaks = c('less_steep', 'as_steep'), values = apal[c(4,5)], 
                         labels = c('Less steep', 'As steep'),
                         name = 'Region, as compared \nto maximum derivative, \nis...') + 
      jftheme)

    cat(paste0('\n\n### Smooth estimate\n\n'))
    
    mf_points <- NULL
    if(include_points){
      mf_points <- geom_point(data = mf, aes_string(y = seven_nets[[i]]), alpha = .15)
    }
    
    #This is where we finally plot the smooth! So it's the correct place to add
    #in the mean of the response. We have to add it to both the line and the ribbon.
    print(ggplot(asmooth[aderivative[, .(Age, compared_to_max)], on = 'Age'],
           aes(x = Age, y = est + response_mf_mean)) + 
            mf_points +
            geom_ribbon(aes(ymin = lower + response_mf_mean, ymax = upper + response_mf_mean), alpha = 1, fill = apal[2]) +  
            geom_line() + 
            geom_line(aes(color = compared_to_max, group = 1), size = 2, alpha = .8) + 
            scale_color_manual(breaks = c('less_steep', 'as_steep'), values = apal[c(4,5)], 
                               labels = c('Less steep', 'As steep'),
                               name = 'Region, as compared \nto maximum derivative, \nis...') + 
            labs(y = 'est') + 
            jftheme)
}

mean_wholebrain_T1wT2w

Derivative

Age where the median posterior derivative is at its max is marked.

Age at maximum derivative

Difference from derivative at steepest age

Smooth estimate

Yeo1_myelin

Derivative

Age where the median posterior derivative is at its max is marked.

Age at maximum derivative

Difference from derivative at steepest age

Smooth estimate

Yeo2_myelin

Derivative

Age where the median posterior derivative is at its max is marked.

Age at maximum derivative

Difference from derivative at steepest age

Smooth estimate

Yeo3_myelin

Derivative

Age where the median posterior derivative is at its max is marked.

Age at maximum derivative

Difference from derivative at steepest age

Smooth estimate

Yeo4_myelin

Derivative

Age where the median posterior derivative is at its max is marked.

Age at maximum derivative

Difference from derivative at steepest age

Smooth estimate

Yeo5_myelin

Derivative

Age where the median posterior derivative is at its max is marked.

Age at maximum derivative

Difference from derivative at steepest age

Smooth estimate

Yeo6_myelin

Derivative

Age where the median posterior derivative is at its max is marked.

Age at maximum derivative

Difference from derivative at steepest age

Smooth estimate

Yeo7_myelin

Derivative

Age where the median posterior derivative is at its max is marked.

Age at maximum derivative

Difference from derivative at steepest age

Smooth estimate